/********************************************************************
Name: Table A10: Effect of Dem Sheriff on Pre-Treatment Detainer
	Compliance Rate
Author: Dan Thompson
Date: August 2019
Notes:
********************************************************************/

* Set the working directory 
gl root = "$sanctuary_root"
gl root = "YOUR_WORKING_DIRECTORY"


/*
Simulate treatment effects under two different simple models
*/

* Bring in the RD analysis dataset
use "$root/data/rd analysis data.dta", clear
keep if policy_year==2018
set seed 73890412

* Set up a container for the simulation output
tempname output
postfile `output' sim_type sim model bw b se n ///
	using "$root/data/trump era rd policy simulation results.dta", replace

* Define a list of functions of the running variable, controls
local func_rv1 = "treat_rv"
local func_rv2 = "rv2 rv3"
local func_rv3 = "rv2 rv3 treat_rv treat_rv2 treat_rv3"
local func_rv4 = "rv2 rv3 rv4 rv5"

* Set up the bandwidths I'll use across models
local bw1 = 10
local bw2 = 25
local bw3 = 25
local bw4 = 25

* Set up the simulation variables
gen sim_enforcement = .
gen rand = .

* Find the number of sheriffs in 287g and listed as sanctuaries
count if enf==1
local enf = r(N)
count if enf==-1
local no_enf = r(N)

* Run all of the regressions
qui forval i=1/1000 {
	forval s=0/2 {

		* Take a draw from the hypothesized distribution 
		replace sim_enforcement = 0
		replace rand = runiform()
		sort rand // assuming no effect, anyone can join 287g
		if `s' gsort treat rand // assuming separation, only reps can join 287g
		replace sim_enforcement = 1 if _n<=`enf' // add the flags for the random subset
		if `s'==2 replace sim_enforcement = 1 if _n<=2*`enf' // add the flags for 2x the random subset
		gsort -rand // assuming no effect, anywhere can be a sanctuary 
		if `s' gsort -treat -rand // assuming separation, only dem can provide sanctuary
		replace sim_enforcement = -1 if _n<=`no_enf' // add the flags for the random subset
		if `s'==2 replace sim_enforcement = -1 if _n<=2*`no_enf' // add the flags for 2x the random subset
		
		* Loop over the five estimators
		forval m=1/5 { 
			
			* Run the OLS regressions
			if `m'!=5 {
				reg sim_enforcement treat rv `func_rv`m'' ///
					if abs(rv)<`bw`m'', vce(clust election_id)
				post `output' (`s') (`i') (`m') (`bw`m'') (_b[treat]) (_se[treat]) (e(N))
			}
			
			* Run the rdrobust routine
			if `m'==5 {
				rdrobust sim_enforcement rv, vce(cluster election_id)
				post `output' (`s') (`i') (`m') (`e(h_l)') (e(tau_cl)) (e(se_tau_cl)) (e(N_h_l) + e(N_h_r))
			}
		}
	}
	noi di `i'
}

* Hang onto the regression results
postclose `output'

* Create a summary file based on the simulation results
use "$root/data/trump era rd policy simulation results.dta", clear


/*
Run the regressions to get the estimated effects
*/

* Bring in the analysis data
use "$root/data/rd analysis data.dta", clear
keep county_id election_id year pre_year policy_year ///
	share_detained_sheriff pre_share_detained_sheriff ///
	all_sheriff_per1k_pop pre_all_sheriff_per1k_pop ///
	enforcement in_287g sanctuary ///
	treat rv rv2 rv3 rv4 rv5 treat_rv treat_rv2 treat_rv3

* Set up a storage variable for the output
gen est_b = .

* Run the regressions and store the results
reg enforcement treat rv treat_rv if abs(rv)<10
replace est_b = _b[treat] if _n==1
reg enforcement treat rv rv2 rv3 if abs(rv)<25
replace est_b = _b[treat] if _n==2
reg enforcement treat rv rv2 rv3 treat_rv treat_rv2 treat_rv3 if abs(rv)<25
replace est_b = _b[treat] if _n==3
reg enforcement treat rv rv2 rv3 rv4 rv5 if abs(rv)<25
replace est_b = _b[treat] if _n==4
rdrobust enforcement rv
replace est_b = e(tau_cl) if _n==5

* Save the five estimates
keep est_b
keep if est_b!=.
gen model = _n
tempfile real_estimates
save `real_estimates'


/*
Calculate the probability that either model is the true model
*/

* Bring in the simulated and real estimates
use "$root/data/trump era rd policy simulation results.dta", clear
keep model sim_type sim b
merge m:1 model using `real_estimates', keep(1 3) nogen

* Set up the containers for the density estimates
forval s=0/2 {
	gen density`s' = .
	gen est_density`s' = .
}

* Get the density of effects under different models by regression type
qui forval m=1/5 {
	forval s=0/2 {	
		sum b if sim_type==`s' & model==`m'
		replace density`s' = normalden(b, r(mean), r(sd)) if model==`m'
		replace est_density`s' = normalden(est_b, r(mean), r(sd)) if model==`m'
	}
}

* Get the distribution of likelihood ratios under the two nulls
* and get the likelihood ratio for the estimated effect under both nulls
forval s=1/2 {
	gen LR`s' = density`s'/density0
	gen lr`s' = log(density`s') - log(density0)
	gen est_LR`s' = est_density`s'/est_density0
	gen est_lr`s' = log(est_density`s') - log(est_density0)
	gen prob_z`s' = 1/((1/est_LR`s') + 1)
	gen prob_z0_`s' = 1/(est_LR`s' + 1)
}

* Make a dataset of the p values
keep sim sim_type model est_b LR1 lr1 est_LR1 prob_z1 prob_z0_1 ///
	est_lr1 LR2 lr2 est_LR2 est_lr2 prob_z2 prob_z0_2
reshape long lr est_lr, i(sim sim_type model) j(sim_type_lr)
keep if sim_type_lr==sim_type
gen p = lr<est_lr
collapse p prob_z1 prob_z0_1 prob_z2 prob_z0_2 est_b, by(sim_type model)

* Print out the table
list model est_b prob_z0_1 if sim_type==1

